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FOREWORD 


Applied  math  uaticians  are  ocrti  nu-.-usIy  concerned  with  numerical 
methods  to  obtain  answers  t.o  pi- at  tioal  problems.  Many  of  these  numerical 
methods  revolve  arouna  a lull  hv»r  foimule  for  polynomial  representation 
of  a function  due  to  Lagrange  • It  is  the  purpose  of  this  paper  to  present 
this  formula  and  apply  it  to  many  numerical  methods  and  to  collect 
together  appropriate  tables  which  will  reduce  the  amount  of  labor 
necessary  to  obtain  numerical  answers 

Although  this  formula  should  only  be  applied  to  functions  which  may 
be  approximated  by  a polynomial  it  is  recalled  that  Weierstrass  Theorems 
shows  this  property  to  hold  *‘c  r a large  class  of  functions  It  is,  of 
course,  true  that  the  degree  of  the  polynomial  is  not  specified  and  if 
a function  is  fitted  by  a polynomial  to  a certain  specified  degree,  the 
accuracy  of  the  fit  should  be  checked  . 

The  collection  of  derivation  of  these  formulae  and  tables  form  a 
part  of  a study  of  numerical  methods  and  the  work  was  performed  under 
Allotment- 45007 , Job  Ordor  2.G  -13  . 

K.  L-  Niel_en 

Head,  Mathematics  Division 
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i INTRODUCTION 


One  of  the  most  applicable  formulas  in  numerical  analysis  is  one 
which  passes  an  degree  polynomial  through  (n  f 1 ) given  points. 

This  formula  is  credited  to  Lagrange  axid  is  usually  referred  to  as 
Lagrange’s  interpolation  formula  sfee  ts  prime, rv  use  was  for  inter- 
polation problems.  However,  it  car,  be  used  as  the  basis  for  many  other 
numerioal  analysis  problems  and  in  special  cases  tables  of  coefficients 
may  be  oalculated  which  reiuce  the  problem  to  simply  one  of  obtaining 
the  sums  of  products  It  is  the  purpose  of  this  paper  to  present  the 
complete  development  of  the  formulae  and  list  some  special  tables. 

2 THE  FUNDAMENTAL  FORMULA 

Let  there  be  given  the  values  of  the  ordinates  y , y^  , . ..  , y^ 

of  the  function  y = f(.i)  at  fie  fr:  -*•  1 ) points  

The  polynomial  of  the  n—  degree  through  these  points  can  be  written 
in  the  form 


SI-53-35 


If  we  let 


(2)  P(x)  = (x  - xq)(x  - x x ) 


'x  - x ) = n (x  - x . ) 

n j=0  J 


p (x)  = (x  - x.  y1  p(x)  = (x  - x y1  n (x  - x ) 
11  1 j=0  J 


then  formula  (1)  becomes 


P (x)  Px(x) 


L(I)  = ^ + *1 


P (x) 

n 

P (x  ) yn 
n'  n 


It  is  easily  seen  chat 


P (x.  ) = 0 if  i £ r 
r ' l 


so  that 


P,(*  ) 


r r / \ 

■/  ■ I y„  = y {r  = o , . . . ,n) 


L(lr>  - FTT7  Jr  “ yr 
r r' 


and  the  polynomial  given  by  (1)  is  one  which  passes  through  the  (n  + 1) 
points  (xr,yr),  (r  = o,...,n)  In  general,  P^(x)  is  a polynomial  of 

degree  n ; i .e . , 


_ / v n r 

P. (x)  = a.  x + a.  x 
iv  ' i,n  i,n-i 


+ . . + a.  with  a.  =1 

i,o  m 


_ / v n n-i 

P.  (x,  ) = a.  x.  + a.  x.  + . 
i ' i ' i ,n  i i ,n~i  i 


+ a.  = k. 

1,0  1 
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It  is  desired  that  L(x)  approximate  tne  function  f(x);  i-e., 

(8)  f(x)  - L(x)  + R(x) 
where  the  remainder  R(x)  is  such  that 

(9)  R{x^  ) = 0 for  (i  - 0,1,  . .,r») 

Furthermore , if  we  let 


R'O  = PU) 


Q(x 


\ 

) 


and  consider  any  function 

(10)  cp(z)  = £(z)  ••  L(z ) - P(z)  Q(z) 

such  that 


(p(x^)  = 0 and  cp(t)  = 0 

if  x ^ x^  for  all  i = 0,1,. -,,n;  then  <p(z)  vanishes  at  n + 2 

points.  By  repeated  application  of  Rolle's  theorem  <pn+1(£)  vanishes 

where  x < £ < x . However, 
o — - n 

Ln+1(z)  = 0 aDd  Pn+1(z)  = (n  + 1)'. 

so  that  upon  differentiating  equation  (10)  we  have 

0 = fI1+1(C)  - (n  + 1)1  Q (x ) 
or 
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/ 


Q(x) 


1 

(n"+T)  1 


,n-r  1 
f 


U) 


and  at  z = x equation  (10)  yields 


(11)  f(x)  = L(x)  + P(x)  fn^(C)/(n  4 l)< 


Thus  we  may  write 


(12) 


f(x)  = .1 
1=0 


P,(x) 


P(x) 

rn  + 1 ) ' 


fn4l(?) 


where  £ lies  in  the  interval  [x  ...  x ). 
v on 

The  last  term  of  formula  (12)  is  in  a sense  a remainder  term  and 
is  a measure  of  the  accuracy  of  the  fit  of  the  polynomial  L(x) 

Formula  (4)  is  invariant  under  a linear  transformation.  Let  us 
make  the  transformation 


x = hu+a  x . = h u . + a 

l i 


then 


P.(x)  = (X-X.)"1  P(x) 


(x-xi)1  (x-xo)(x-x1)  (x_xn) 


= (hu  + a -hu^  - a ) (hu  + a - hu^  - a ) (hu  + a - hux  - a ) . . . . (hu  + a - hun~a ) 

=h  1(u  - u.)  1 hn(u-  u )(u  - u.  ) ....  (u  - u ) 

\ x ' v o 1 n' 

= h^~ 1 (u  - u.r1  P(u) 

= h0'1  Pi(u) 
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P.  (x4  ) = (x  - x )(x.  - X,  ) • * - (x.  - x.  )(x.  - X.  ) ...  (x.  - X ) 

i'  i7  ' i o 1 - i i-i  i i+i'  ' i n7 

= (hu.+a-hu  - a)  . . . (hu.+a  hu.  - a)  (hu  + a-hu.  - a ) . . . (hu.+a-hu  - a) 

' i o i i-i  7*  l i+i  7 ' l n 7 


, n- 1 , w 

h (u.  - u ;>u  “ u 
' 1 0*1 


(u.  - u.  }(u  - 0.  ) • 

• 1 1 -l7'  1 1+1 7 


(u.  - u ) 

1 XI 


, n-i  _ , x 
h P.  (u.  ) 

l l 


Therefore 


n P.  (x  ) 
ilo  P- yi 


R P.(u) 

i 5—7 — v y = L,(u^ 

i=o  Pi(u.)  '’l 


3.  EQUALLY  SPACED  INTERVALS 


If  the  values  of  x.  (i  = 0....,n)  are  giver  at  equally  spaced 
intervals , we  have 


x = x + i h 
1 o 


where 


h = A * = x , - x = x - x , f ~ = ],...  , n ) 
i o r r- 1 


and  also  let 


x - x 


or  x = x + uh 
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The  polynomial  P^(x)  now  takes  a special  form 


_!  n 

(15)  P (x)  = (x  - x ) .n  (x  - X ) 

1 1 J=0  J 


1 


= (x  + uh  - x - ih)  FI  (x  + uh  - x - jh) 


J — o 


h 1 (u  - i ) 1 n h(u  j ) 

J — o 


= hn(u  - i)  n (u  - j ) 

,i=o 


Now 


(16)  n (u  - j)  = u(u  - 1)  • 

J«o 


(u  - n) 


where 


n+i  n+l  _n+i  n 
S u + S,  u + 
o 1 


„n+l 

+ S u 
n 


_n+l  / . 
si  ^ 


= 0,...,n)  are  Stirling's  Numbers  of  the  first  kind. 


These  numbers  have  the  property  that 


(17) 


n+l  n n 

S.  = S.  - n S.  , 

l l l-l 


and  this  recurrence  formula  is  easily  used  to  build  up  a table  of  these 
numbers  (see  Table  I). 

For  equally  spaoed  intervals  we  have 


6 


(18) 


P.  (x. ) = 
1 1 


n 

n 

j=o 


! x . . - -X  ; . . 
i J J+i 


= h‘ 


n 

n 

j=o 


(i  - j ) 


so  that 
(19) 


and 


= hn  {-if  1 


Pt(*) 

P.’(x  ') 

L * 1 


h“  n (u-j) 


j^o 


b it  (n-  i }' 


•>**  . . >o-i 

(-0 


= ( - 1 )n  1 [ i ( n - i)‘ ] 1 n (u  - j) 

J=o 


(20)  L(x)  = Lr  (u)y  * L (u)y1  + . . + Ln(u)yn 

where  L^(u)  is  given  by  (19) 


Since  L^(u)  are  functions  of  n and  u it  is  possible  to  compute 

tables  of  these  coefficients. and  thus  the  value  of  L(x)  for  a 
specified  x can  be  computed  by  a sum  cf  products  In  order  to  cut 
down  the  range  of  values  on  u for  such  a table,  a shift  to  a "center" 
value  is  made  which  changes  the  notation  somewhat  Let  us  remember 
the  n + 1 point  x , x2  . . . , in  the  following  manner. 

Case  1.  Even  number  cf  points:  n + 1 = 2r 

In  this  case  we  could  remember  the  points 


-r 


-r+i 


x 

r 


and  we  are  faced  with  a choice  of  picking  the  "center"  value  to  be 
either  x or  x . The  general  practice  is  to  pick  x i to  be  the 

"center"  point  and  call  it  xq  The  numbering  then  becomes 
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X ,X  , X , X , X_,  . . . , X 

-r+i  -r+a  -i  o i a r 


Example • Given  6 points,  they  are  numbered 

x~a ’ x-i’  V V Xa ’ x 3 

We  now  have 


x.  = x + ih 
1 0 

(i  - -r+i , . . 

and 

x = x + ph 
0 

so  that 

(21) 

P.(x)  = (x  - X.  )' 

n 

1 n (x.  - x . ) 
j=0  J 

= h 1 (p-i 

r 

)” 3 n (x  + ph  - X - j 

j=-r+i  0 0 

= h2r‘:  (p  - i)“'“  n (p  - j)  . 

j=-r+i 


The  nroduct 


n (p-j) 

j=-r+l 


= (p+r-l ) (p+r-a ) . . (p+i ) p (p-i ) . . . (p-r ) . 


n r r 

(22)  P.  (x.  ) = n (x.-x  ).  , = h r 3 n (x.  +ih-x  -jh)=h  r 1 II  (i-j).,. 

11  j=o  1 J J*1  j=r+i  0 0 j-r+i  1+J 


= h2r  l [ (j+r-1 ) (i+r-a  ) . . . (i+i  ) (i  ) . . . (i ) (-1  )(-a  ) * - - (r— i ) } 


h2r  1 [ (i+r-i  )'  ( r-  i )?.][- 1 ] r 1 
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and  the  coefficients  of  y become 


(23)  A.  (p  ) = (-1  )*  A [(i  + r-l)*{r-i)t  ]_1  n (p  - j ) . . 
1 j=-rn 


In  the  tabulation  of  these  coefficients  the  following  property  is  useful 


(24) 


Ai(p)  = (1-  P) 


Case  2.  Odd  number  of  points;  n t 1 = 2r  + 1 
In  this  case  the  poirts  are  numbered 


X ,X  , . . . , X . X , X , . . . , X 

-r  -r+i  -loi  r 


and 

(25) 


A.  (p ) = (-1 ) [ (i  + r)'.  (r-  i )■  ] n (p-j).  • 

i i=-r  J-1 


with  the  property 


(26) 


Af  (p ) = A ^(-p) 


Let  us  consider  as  an  example  the  case  of  5 points,  and  let  p = .7. 

It  is  required  to  find  A , A , A , A , A . Substituting  into 

-a  -i  o l a 

formula  (25)  we  obtain 

A_a(.7)  = (-l)2+a[(-2+2)*(2+2)']'1(.7  + l)(.7+0)(.7-  1)(.7-  2) 

= (-1)4[(4)*]"1  [(1  7)(.7)(-.3)(-1.3)] 

-si  < 4641» 

= .0193375 

9 


W-53-35 


A ^.7)  = (-l)a+1[(-l+2),-(2  + l),r'1i  ”,  7 +2)(.7  -»  0)(.7  - 1)(.7  - 2)] 
= - (3*.  )"  1[  (2.7  ) ( 7 )(-  3 ) - 1 3 ) ] 

= -.12233. 


A0(.7)  = (-l)a[(2!)(2l)]'1[(  74  2)(.7  + l)(  7-l)(.7-2)J 

= .447525- 

A1(.7)  = ( — 1)1[  (3*  )(1*.  )]_1(  7 + 2)(.7  + 1 ) ( . 7 + 0)(  7-2) 

= .69615- 


Aa(  7) 


(-1 )°[4*  (O'.  ) ] ( . 7 + 2 ) ( 7+l)(.7+0)(.7-l) 

- .0401625. 


Extensive  tables  of  these  coefficients  have  been  published  (see  ref- 
erence (a)).  A limited  table  for  five  points  is  given  in  Table  II 


4.  INTERPOLATION 


One  of  the  prime  applications  cf  Lagrange's  formula  is  in  interpo- 
lation. The  problem  is  to  find,  a value  of  a function  at  a point  x 

which  falls  between  tabulated  values  at  points  r.  , x . - . . 

ci 

A.  EQUALLY  SPACED  INTERVALS 

The  simplest  case  is  the  one  in  which  the  tabulated,  values  are 
given  at  equally  spaced  intervals,  then  formula  (20)  gives  the  desired 
results.  Furthermore,  the  coefficients  may  be  tabulated  and  the  problem 
is  reduced  to  simply  one  of  finding  the  sum  of  products  of  two  numbers 
and  is  easily  performed  as  one  operation  on  a calculating  machine.  For 
this  purpose  the  coefficients  are  tabulated  about  a "midpoint"  and  we 
have 
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(27) 


t(*)  = A.ry.r 


+ A 


•o  -rti 


+ A 


7 -*■  A.y. 

c o 11 


. + A y 
r r 


The  procedure  is  explained  by  an  example. 

Example  1. 

Find  f(1.77)  by  a 5-poir.t,  Lagrangi  ar.  formula  from  the  table  of 
values 


xi 

1.5 

1.6 

1.7 

1.8 

19  2.0 

yi 

48  .09375 

65 . 53600 

87,69705 

115-47360 

149.86915  192.00000 

A i ( - 7 ) 

.C193375 

- .1228500 

4475250 

.6961500 

- C401625 

A,  (-3) 

- .0261625 

.2541500 

..8895250 

-.1368500  .0193375 

o 1.77-1.7 

Solution  Let  x =17.  then  p = — = = = .7 

o h .1 

The  Lagrangian  coefficients  hi")  are  obtained  from  Table  II 
and  listed  above. 


f (1  77) 


2 

2 A. ( 7)y  = 105.49336 

■i  _ i l 

1=~2 


If  we  had  let  x =1,8, 
o 


1 .'77  - 1.8 
p = , 


we  obtain 


2 

f(l.77)  = 2 A (-.3)y.  =106.49348 

i=-2  1 1 


B.  UNEQUALLY  SPACED  INTERVALS 

The  case  of  unequally  spaced  intervals  makes  the  method  more 
difficult  to  handle  tut  increases  its  importance  since  fewer  methods 
are  now  available.  It  now  becomes  necessary  to  return  to  formula  (4). 
A schematic  may  be  devised  which  greatly  aids  the  computation.  First, 
a linear  transformation  is  made  on  the  given  . (i  = 0....,n)  in 

order  to  obtain  small  integers  in  so  far  as  is  possible. 
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Form  now  the  square  array 


r 


f 

x - 

X 

0 

X 

0 

xi 

X 

c 

- xa 

...  X - X 

c n 

j 

xi 

- X 

0 

x - 

xi 

X1 

- xa 

X , - X 

1 n 

(28)  J 

xa 

- X 

0 

X2 

- xi 

X - 

xa 

X_  - X 
3 n 

X 

n 

- X 

0 

X 

n 

- xi 

X 

n 

xa 

...  x-x 

n 

V 

We  note  that  the  product  of  the  principal  diagonal  is 


(29)  P(x)  = (x  - xo)(x  - xx)  ■ (x  - xn)  . 

The  products  of  the  elements  in  each  row  yield 

(30)  R.  = (x.-x  ) ...  {x-x  ) ...  (x.  - x ) = (x-  x.  ) P.  (x.  ) 

l i o l ' l n * i'll 


Thus 


(31) 


P.  (x)  (x  - x. ) Xp(x)  , , 

1 ' 1 ' ‘ p.x) 


p“Tx~7  ~ 7 

1 ' 1 ' (x  - X . ) 

■J  * 


R. 


and 


(32) 


P(x)  P(x) 

+ —5  ~ y3  + • • • ~ 5 ~ 
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The  square  array  (28)  may  now  bo  augmented  by  x-wc  columns  R^  and 

y i / Hi  (for  speedy  operation  R^  need  net  be  recorded)  The  sum  of  the 

last  column  multiplied  by  the  product  of  the  terms  of  the  main  diagonal 
yields  the  desired  result. 

Example  2 . 


Given  the  values 


X 

0 

3 

9 

12 

15 

21 

27 

y 

150 

108 

0 

- 54 

100 

144 

-34 

Find  y at  x = 18  - 


Solution 


We  first  make  the  transformation  x = 3s,  then 


s . 

0 

~T 

3 

4 

5 

7 

9 

l 

and  s = 6.  The  array  of 


. - s.  becomes  after  augmented  by  R.  and  y /r 


s . and 
i 


s . -s 
1 0 

VS1 

Vs* 

s . -s„ 
1 3 

VS4 

s.  -sc 
1 6 

s . -sc 

1 6 

R. 

l 

6 

-1 

-3 

-4 

•5 

- 7 

- 9 

22680 

.006613757 

1 

5 

-2 

-3 

-4 

-6 

-8 

-5760 

- -01875 

3 

2 

3 

-1 

-2 

-4 

6 

864 

0 

4 

3 

1 

2 

-i 

-3 

- $ 

-360 

,150 

5 

4 

2 

1 

i 

-2 

-4 

320 

-.3125 

7 

6 

4 

3 

2 

- 1 

-2 

2016 

- .071428571 

9 

8 

6 

5 

4 

2 

-3 

-51840 

.001620370 

The  product  of  terms 

in  principal 

diagonal  P(s 

)=  540. 

- 244444444 

y = (540)(  .-244444444)  = -132. 
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G INVERSE  INTERPOLATION 


Lagrange's  formula  adapts  itself  very  nicely  to  inverse  interpolation 
since  formula  (2)  is  simply  a relation  between  two  variables,  either  of 
which  may  be  considered  the  independent  variable.  Thus,  we  can  write  x 
as  a function  of  y 


(33) 


L(y)  = 


V*) 

+ rm 

n n 


x 

n 


Normally,  the  values  of  y will  be  unequally  spaced  and  the  mothod 

just  described  must  in  general  be  employed  Furthermore,  it  will  be 
more  difficult  to  find  a linear  Irans forma  tier  which  will  reduce  the 
size  of  the  numbers.  Let  us  consider  an  example 

Example  3. 


From  a five-place  table  of  natural  sines  we  have 


X 

30 

31 

33 

34 

36 

y 

.50000 

,51504 

54464 

55919 

7','  9 

Find  x when  y = .52992- 
Solution 


Let  y = .Is  , then 


sL  5.0000  5.1504  5 4464  5.5919  5-8779  5.2992 


and  the  computational  form  is 


s.-s 
i o 

Vsi 

Vs* 

s- 

1 3 

VS4 

R 

i 

x./R. 
i l 

.2992 

-.1504 

- .4464 

- 5919 

8779 

.010438234 

2874.049 

.1504 

.1488 

- .2960 

-.4415 

- .7275 

- 002127679 

-14569.867 

.4464 

.2960 

-.1472 

-.1455 

-.4315 

- .001221146 

-27023.796 

.5919 

.4415 

.1455 

-.2927 

-.2860 

.003182957 

10681.891 

.8779 

.7275 

.4315 

.2860 

- .5787 

- .045611921 

-789.267 

Sum  = -28826.990 


P(x)  = - 001110065;  x = 31  9998. 
To  five  place  accuracy  x = 32 . 
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D . MULTIPLE  INTERPOLATION 

Interpolation  in  tables  of  multiple  arguments  is  a laborious  pro- 
cedure and  is  usually  accomplished  by  repeated  application  of  interpo- 
lation on  a single  variable  Repeated  application  of  Lagrange's  Formula 
can  best  be  illustrated  by  example-  VJe  shall  limit  ■•■ur  attention  to 
equally  spaced  intervals  , for  unequal-)  v spaced  intervals  repeated  applica- 
tion of  the  technique  of  4 B must  be  employed. 

Consider  y = f(x,r)  and  suppose  the  table  of  values  is  arranged 
as  follows 


< X < X. 


Let  it  be  desired  to  find  a value  of  y(x  r)  where  a2  - A - *.3 
rx  < r < r2  by  a three  point  Lagrangian  formula 


r - r. 


We  obtain  p = 


and  p = and  utilize  the  points 

Ax  r Ar 

x,  , , x,  and  r , r,  , r_ , denoting  the  Lagrangian  coefficient  by 

X 4 J Q a 4 

xxx  r r r 

A , , A , A and  A , A , A.  If  we  interpolate  first  with  respect 
— A O 1 ~ X 0 -a. 

to  x at  each 


r.  (i  — 0.1. we  have 

l 
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(34) 


X 

< 

II 

y 

-f 

AX 

y 

xo 

___  1 

10 

0 

20 

xl 

II 

> 

X 

y 

+ 

kX 

y 

-1 

'n 

0 

' 21 

= A* 

y 

+ 

•J 

A 

y 

X2 

-l 

' 1 2 

0 

7 22 

Y 

A y 


32 


Then  interpolating  with  respect  10  r yields 


(35)  y(x,r)  = A , y + A y , * A'  y 

. r , . x , x x , 

= A A y - A v A y J 
-l  1 -i  'lo  o '20  l ;30 

„ r r „ x x.  , x i 

+ A [A  v -*  A v + A V J 

0 -X  11  O 21  1 31 

. r r . x . x x , 

+ A [A  y t A v f A y J 

1 1 ' 1 2 O '22  1 32 


nr  AX  „r  „X  „X 

A Ay  y A Ay  + A Ay 
-1  1 710  -1  0 20  -1  1 


30 


r . x 


r . x 


+ A A y +AAy.  +AAy 
o -i  'n  o o 'ai  o i '31 


r x r x r*  Y 

+ A A y - A A y * A'  A y 

1 -1  13  1 0 22  1 1 ' 


32 


It  is  seen  that  nothing  is  gained  by  attempting  to  use  the  last  expres- 
sion since  this  would  require  the  recording  ol’  the  9 products 

j 

A.  A.  (i,j  = -1,0,1)  while  using  the  first  expression  only  4 recordings 
i J 

are  necessary.  The  second  method  would  have  application  only  in  the 
special  case  where  p^  and  p^  remain  constant  for  a large  number  of 

interpolations.  The  most  expeditious  method , therefore , is  to  compute 
the  three  y (i  = 0,1,2)  and  then 


(36) 


p r 

y(x.r)  = A y + A y «-  A y 
' -l  x.o  c XL  i x 2 
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5.  DIFFERENTIATION 


Lagrange's  formula  may  be  used  to  tind  the  derivative  of  a function 

only  at  discrete  values  x fi  = 0,  n) 

t ' f)  = L.S)  + R ( x ) 

= f (x)  = L (x)  + R (x. ) 

= P<r)  f(n+l)(C)/(n  4-  1)! 

(38)  R'(x)  = P"(x)  f(nJ’‘)(C)/(n  4-  l)'  + p(x)  ffn+a)(?)/(n  + l)i 

The  second  term  of  this  expression  would  be  difficult  to  find  even 

/ t 

if  it  is  known  that  f''h+2^(x)  exists  However,  if  it  were  evaluated 
at  a given  point,  x^,  then  P(x^)  - 0 and  only  the  first  term  remains. 

In  practice  the  remainder  term  is  only  used  to  check  the  accuracy. 

The  derivative  of  L"(x)  is  given  by 


y = f(x)  which  is  known 
From  formula  (8)  we  have 


so  that 
(37) 


dy 

dx 


Since  by  formula  (11) 

R (x ) 


(39) 

and 


L (x ) = L (x)y  + L'(x)v  + . . + L'(x)y 

' ' o o i • ’ l n n 


(40) 


L'(x) 


1 d r , » 
P.  (x.  ) dx  Pi^ 

l i 


PTF7 

ii 


[ (x-x 


(x-x.  )(x-X.  , ) 

v 1-1  1+1 


X-X 


x-x.  X-X. 

1-1  1+1 


+ . . . + 


X-X 


W-53-35 


(40) 


K<‘>  = vrr 


n 

1 

[2  1 I 

n (x  - x . ) 

j=o  J 

i-o  x x 
J-°  j _ 

J41 


Pi(x) 

FTT7 

i i 


2 -i- 

j=o  ' \, 


j+i 


Formula  (39)  gives  an  expression  for  finding  the  derivative  of  a function 
at  a general  value  of  x.  The  computation  of  the  coefficients  of 

(i  — 0 n);  namely,  L/(x).  as  given  by  formula  (40)  is  somewhat 

difficult  to  manage.  It  can  best  be  handled  by  lorming  the  square  array 
(28)  and  computing  P (x)/P  (x^)  as  was  done  by  formula  (31). 


A.  THE  DERIVATIVE  AT  xk  FOR  UNEQUALLY  SPACED  INTERVALS 

Of  special  interest  is  the  derivative  evaluated  at  one  of  the  given 
points  x (i  = 0 , ....  n).  Formula  (40)  takes  on  two  forms: 

n 

(a)  k 4=  i • In  this  case  the  product  H 'x~  x.).,.  = 0 except 
j_0  J JT1 

when  it  is  multiplied  by  (x^  - x ) 1 thus 


(41) 


L'(x) 


n 

n 

j=o 


jkk  J 


- — 7 r [ (x,  -x  ) 

P. (x. ) 1 ' k o' 
i l 


=— r [ (x,  -X  ) - (x.  -X.  Vx.-X.  ),.  (x, -x  )]/(x. -x.) 

P.  (x.  ) k o'  x k k-l  k k+r  k n'  v k i' 


l l 


p7^-7  (xk  - Vo*  (xk  - V1 
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If  we  let 

(42) 

we  have 

(43) 

(b) 

(44) 
and 

(45) 

Thus  the 

(46)  I 


D..  = (x.  x.  ) P. (x.  ) , i * k 

lie  V l l l > ' 


Li<*>  - 4 & - Vj*  • 


k = i . In  this  case  wo  have 


n (x.  - x . ) . . = P (x.  ) 
i J Jfi  i i 


P (x.  ) n , 

, ■ > \ i i - 1 

L.  IX  1 = — , 2 

1 P. (x.  J ,_r  x.  - X . 

1 ]•  J=°  1 J 

jt=l 


= 2 (x.  - x . ) 

complete  expression  fox'  l/vx^)  becomes 

/ (x  ) = l/(x  )y  + ...  + LT(x  )y  + ...  + L'(x  )y 
k o k o k k k n k n 


y0  n yx  n 

p n ’VX.,  - x.).,,  + n (x,  - X.).  . + .. 

D0k  j=°  k J <*k  Dik  j~o  k J Jfk 


n 


+ y.  2 (x  ,-x  . ) * 
k j=o  k j'j*k 


yn  n 

• • + rf“  n (xk-xi 

Dnk  J=0  k J J*k 


( vxj  >j+k 


Dok  Dik 


nk 


+ Yk  jS0  ^Zk'XPj+k 


19 


W-53-35 


or 


n 

(47)  1/  (x)  = n 

j=o 


< V*j  >0** 


n 

Z 

i=o 


'ik 


l+k 


4 j!0 


This  formula  for  the  derivative  of  y = f':<)  at  i = t can  be  found 

by  the  use  of  a schematic  similar  to  the  one  given  in  Section  4.B.  Form 
the  square  array 


X 

X 

X 

x_ 

X 

x, 

X 

X 

0 

1 

0 

2 

0 

k 

0 

n 

X 

X 

X 

X 

...  X. 

X. 

...  X 

X 

1 

0 

1 

2 

k 

1 

n 

X - 

X 

X 

X 

. . X 

...  X 

X 

2 

0 

2 

1 

2 

k 

2 

n 

X,  - 

X 

x,  - 

X 

X,  - 

X 

x.  - 

X 

k 

0 

k 

1 

k 

2 

k 

n 

x - 

X 

X - 

X 

X 

X 

X 

n 

0 

n 

1 

n 

2 

n 

k 

Notice  that  no  element  appears  on  the  main  diagonal.  Now  the  product 
of  the  elements  of  each  row  is  P^(.«i).  If  this  be  again  multiplied  by 

(x^  - x^)  we  have  -D^.  Thus  we  can  augment  this  array  by  two  columns, 

D.,  and  y.  D.,  ; the  first  being  obtained  by  multiplying  the  product  of 

l iC  l 1 M ^ ^ 

each  row  by  the  element  in  the  k column  (i.e.,  this  t,j.emen„  twice  in 
the  product)  and  changing  sign  There  will  be  no  entry  in  thb  (i  = k) 

n n 

row.  The  product  n (x  - x.).,,  = n [-(x.  - x,  )].  , is  obtained 

j=o  * J Jt^  j=o  J * Jr* 

th 

by  multiplying  together  the  negative  of  the  elements  in  the  k column. 

The  elements  of  the  second  augmented  column  are  now  summed  and 
multiplied  by  this  product  to  yield  the  first  term  of  the  formula  (47). 


i 


| 
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The  second  term  is  obtained  by  scanning  the  negative  reciprocals  of  the 
terms  in  the  column  as  can  so  easily  seen  from  its  expression  The 

final  value  of  the  derivative  is  then  given  by  formula  (47)  Let  us 
illustrate  the  procedure  with  an  example. 

Example  4. 

Let  us  consider  the  example)  2 given  in  4.B,  and  let  it  be  desired 
to  obtain  the  derivative  at  s.  = 5;  i.e  , x = 12 


s . -s 
1 0 

s . -s 
1 1 

s . -s 
1 2 

s . - s 

t 3 

s . 

1 4 

G "G 

1 5 

3 . - 3 
1 6 

D. 

14 

y.  d'1 

i i4 

-1 

-3 

-4 

-5 

-7 

-9 

18900 

.007936507 

1 

-2 

•7 

— %■> 

-4 

-6 

*8 

-46C8 

- .023437500 

3 

2 

-1 

-2 

•4 

-6 

576 

0 

4 

3 

1 

. ]_ 

-3 

-5 

-180 

.300000000 

5 

4 

2 

1 

■ 2 

-4 

7 

6 

4 

3 

2 

-2 

4032 

- .035714285 

9 

8 

6 

5 

4 

2 

-69120 

.001215278 

n (s 

4 

s.  ) = 320. 
i 

.250 

Z(s4 

-i>‘l 

1 X 

= v + * 

1 

+ 2 + 

i i 

1 " 2 * 

l _ 

4 D 

= 1.2000 

L'(x 

) = (320  ) ( .25  ) 

+ (l.£ 

>)(-100)  = 40. 

B.  THE  DERIVATIVE  AT  xk  FOR  EQUALLY  SPACED  INTERVALS 

If  the  values  of  the  independent  variable  are  given  at  equally 
spaced  intervals,  we  have  from  equation  (20) 

(48)  L(x)  = Lo(u)yQ  + L1(u)yi  + ...  + Ln(u)yR 
where 

(49)  L.  (u)  = (-l)n_i  [i'-(n-  i)')1  U (u-j)...  . 

J=o  Ut1 
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Since  x = uh  + x we  have 
o 


dy  _ dy  du  1 dy 
dx  - du  dx  ~ h du 


and 

(50) 


y'  = l/(x)=i  Z L.  (u)y. 

n i -o  1 *■ 


and 

(51) 


L*r  (u)  = (-!  )n  1 ( i (n  - i ) ' ] ~ n (u-j). 


du 


W 


Furthermore  [see  formula  (40)] 


(W) 


n (u-j) 

J=0 


1 <u  - J >j+i 


J=° 


Hhen  x = x,  we  have 
k 


(53) 


X - X X,  - x 

o k.  o , 

u = -h—  = — — “k 


so  that  at  this  point  formula  (52}  takes  cn  two  special  forms: 

(a)  k =fc  i 

d n 

(54)  n (u-j  ) = lu(u-l).  . (u-k+l)(u-k-l). . . (u-i+l)(u-i-l). . . (u-n)] 


= k(k-l)  1 (-1)  ,.(k-i+l)(k-i-l)...(k-n) 
= (k-i)  : k"  (n-k)!  ( -1 )n  ^ 


22 


and 

(55) 


L'(u)  = (-l)n"1[u(n-  i )'.]'*  (k-  i)"*k!(n-  k)t(-l) 


n-k 


= (-D 


i+k  k '■  ( n - k ) 1 

i'  (n - i )'  (k  - i ) 


Note:  (-l)an'i'k  = (-l)an(-l)-''1+k^  = (-l)1+k 


(b)  k = 


(56)  -r-  .n  (u -i)  = [i(i-l)...(2)(l)(-l)(-2)...(i-n)] 

au  j —O  Jt1 

r1  1 11  1 i 

i 1-1  1 1 i-n 


it  (n  - i (-1  )n_i  2 

0=o 


i - J 


jf1 


so  that 


(57) 


L.'(u) 


= (-D 


n-i  i l ( n - i ) ! ( - 1 ) 
i ' ( n - i ) '• 


n-i 


0=o  1 


j+i 


11 
i 1-1 


1111 
•f  • • 4*  ■ j - «■-=»  — 

2 1 1 2 


1 

n-i 


V 1 ^ 1 

= Z -v  - Z -r 

j=  i J j=i  j 


n-i 


Z y 
j=i+i  0 


if  2i  < n- 1 


i 

Z 

j=n-i+i 


if  2i  > n-1  . 
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Thus  the  coefficients  of  y.  (i  = 0...  ,n)  are  all  functions  of  n,  k, 
and  i,  can  be  computed  once1and  for-  al-i  , and  tabulated.  Furthermore, 
the  table  of  coefficients  is  "negatively  symmetric"  about  the  "midpoint" 
so  that  it  is  necessary  to  tabulate  only  half  of  them.  Let  the 
coefficients  be  denoted  by  , then  for  equally  spaced  intervals 


(56)  y'(xk)  • L-(xk)  = r.  [Ak0  y0  . Ak.  Sj  + 


* A"  y J 

kn  n 


and  Aj\  are  given  by  formulas  (55)  and  (57). 

Now  for  k > in  where  m = ~ n + 1 if  n is  even  and  m = ^-(n  + 1) 

C. 

if  n is  odd,  wo  have 


(59) 


A'  = -A'  , 
ki  n-k,ri-i 


By  formula  (55)  we  have  for  k ^ i 


ki 


= (-D 


i+k 


» f T 


k)’ 


i-  (n  - i)'-  (k  - i ) 


and 


_ / 1x(n-k)+(n-i)  (n-  k)l  (n  [n  k])'. 

n-k  ,n-i  (n  - i )•  (n  - [n  - i 3 ) *•  (n  - k - [n  - ij) 

= (-l)k+i  (n-k)'  (kQ 

(n  - i)'.  (i  ) ••  (i  - k) 


= (-D 


k+i-i 


(n  - k)'.  k*. 

(n  - i )'-  (i  )'.  (k  - i ) 


1 


thus 


-A 


n-k ,n- i 


24 


I 


W 53-35 


» 


For  k = i we  have  by  formula  (37} 


A',.  = Z 

j 


kk  i=ti-  k+l  J 


and 


n--k  ,n-k 


so  that 


A.' , = -A  . , 

kk  n-k  n-k. 


n-n+k  ^ 

2 T 

j=n~k+i  J 


The  tabulation  of  A'^  can  be  made  in  two  ways.  The  first  is  to 

obtain  the  lowest  common  denominator  for  the  A,'.  for  each  n and 

ki 

tabulate  the  coefficients  of  >\  as  integers  with  a final  division  to 
obtain  the  derivative.  Thus 


(60) 


y'(*k>  = 


h D 


n 

Z 

i-o 


where 


D 

n 


and  D is  the  lowest  common  denominator  of  A.'.  ..  Table  III  gives 
n ki 

the  values  of  D and  a'  for  10  points.  The  second  method  • ^o 
n ki 

publish  A'  in  decimal  form,  see  Table  IV 


Example  5 . 

dy 

Find  at  x = 1.7  in  the  data  given  in  Example  1^. 

Here  u = — = 2 = k;  n = 5 h=.l 

• X 

From  Tables  III  (D  =60)  and  IV,  we  he’re 
n 
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yi 

48.09375 

149.86915 

192  00000 

a . 
ai 

3 

-30 

-20 

60 

-15 

2 

A , 
ai 

.050000 

- .500000 

- ,250000 

,033333 

II 

O 

, 

6 

2 a'  y. 
i=o  21  1 

/60  = 248  10650. 

10 


6 

2 

i=o 


yi 


= £48  10658. 


6.  INTEGRATION 


As  in  the  case  of  differentiation  Lagrange's  formula  may  also  be 
used  to  find  the  integral  cf  a function  f(x)  which  is  known  only  at 
discrete  values  (i  = 0,...,n).  That  is  the  integral  of  L(x)  will 

be  a good  approximation  of  the  integral  cf  f{x)  so  that  we  may  write 
b 

(61)  / f (x)dx  = aQyo  + ax  7X  * - + anYn  + R 

a 

where 

1 b 

(62)  ai  = p"T577  f pi'*  ^dx  ( L ~ 0 » • • • i n ) 

i v i ' a 

and 

(63)  R = -fc  P(x)  dx 
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The  problem  of  finding  the  integral  of  f(x)  then  reduces  to  the 
evaluation  of  the  coefficients  a.,  and  the  sum  of  products.  The 

coefficients  a.  depend  upon  the  n + 1 points,  x^  . and  the  limits  of 

integration  (a,b).  If  we  write  p.{x)  as  a polynomial  in  descending 
powers  of  x 1 

(64)  P.  (x)  = + a xr*  1 4 . . . + a.  x + a. 

i i ,n-i  i,i  i,o 


the  integral  may  be  written  in  the  fonn 


(65)  fp.  (x)dx  = b.  xr+1  + b.  x'1  4-  . . 4 b.  r + b. 

l l , n4 1 i ,n  1,1  i,o 


where 


i ,n4-i  n + 1 


(i=0;  . ,n) 


(66)  < b =-?(&  ) (i  = 0 . . ,n  : j = 1 n) 

i , j j ...  j - i 


b^^  = constants  of  integration 


The  coefficients  b. , may 

ij 

by  a row  of  the  constants 


since  they  will  disappear 
Thus 


bo  placed  in  a tabular  form  and  augmented 

P.  (jc.  ) It  is  not  necessary  to  write  b. 
i i ' 10 

when  the  limits  cf  integration  are  inserted. 


n+i 

X 

n 

X 

n-i 

X 

. . . 

x2 

X 

P.  (x.  ) 

l l 

b 

b 

b 

b 

b 

P (x  ) 

o ,n+i 

o ,n 

0 , .U-  1 

0 2 

0 , 1 

0 0 

b 

b 

b 

b „ 

b 

P (x  ) 

n,n4i 

n,n 

r>  ,n-i 

n ,2 

n , l 

n n 

i 
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There  remains  now  to  consider  the  range  of  integration  If  this  range 
changes  so  that  it  is  desired  to  have  uianv  integrals,  the  values  for 
may  be  tabulated  in  the  following  manner 


i 

H 

! x 

i 

o 

M 

x - x_ 
o 2 

V 

*3  Te 

• * * 

0 

ao 

ao 

— 

co 

P 

O 

1 

al 

al 

a 

i 

ai 

■ • • 

2 

a2 

a2 

a2 

• • ' 

a2 

• 

• 

• 

• 

• 

• * • 

n 

an 

an 

an 

For  the  case  of  equally  spaced  intervals 


(67)  P (x)  « P (u)  « u(u-l)  . (u-  i+1 ) (u-i-1 ) . .(u-n) 
and  since  dx  = h du 


n 


2 

j-o 


a. 


i .J 


u 


J 


(68)  /P  (x)  dx  = h/P  (u)  du  . 

The  values  of  a_  may  be  found  easily  by  use  of  Stirling's  numbers  of 
the  first  kind.  Consider  for  example  the  case  of  4 points 

PQ(u)  = (u-l)(u-2)(u-3)  = u3  ~ 6u2  + llu  - 6 ; 

P1  (u)  = u(u  - 2)  (u  - 3 ) = u3  • 5u2  + 6u  ) 

< 


P3  (u)  = u(u-  l)(u  - 3)  = u2 
P3  (u)  = u(u - l)(u  -2)  = u3 


4l?  + 3u  j 
3u2  + 2u 


S 
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If  we  note  that 


P.<u)  = P(u)/(u  i) 

the  polynomial  P(u)  has  for  its  coefficients  the  Stirling  numbers  of 
the  first  kind  and  from  which  we  remove  the  factor  (u  - i).  This  can 
be  done  by  a synthetic  division  and  easily  put  into  a schematic.  The 
removal  of  (u  - o)  does  not  change  the  numbers  so  that  the  first  poly- 
nomial has  the  Stirling  numbers  of  the  first  kind  as  coefficients. 


W-53-35 


W 


are  easily  found  by  the  formula 


P.  (u.)  = i'(n-  i)i  (-l)n_1 

and  realizing  that  P (u  ) = a and  then  they  alternate  in  sign.  A 

o o oo 

column  for  the  lowest  common  denominator  (L.C.D.)  is  attached  in  case 
it  is  desired  to  factor  out  the  denominator;  it  is  to  be  recalled  that 
P^(u^)  is  a denominator. 

Let  us  now  consider  a specific  integral,  say 

/ f(x)  dx 
o 


for  4 points . We  have 


\ = tnW  /p(u) 

o o J 


w[b°‘ 

u + b uJ  + 

b u + b u, 
02  01 

03 

h [i(2>4  - 

! <2>3  + -V1 

(2 )2  - 6 (2) 

= - I [i]  K3>(2>4  - 24<2)3  + 66(2)2  - 72(2>J 

= _ JL  [48  - 192  + 264  - 144 J 
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The  schematic  eases  the  oomputali  or 


For  equally  spaced  intervals  tne  values  of  s . can  be  tabulated 

for  the  range  cf  integration  0-1.  0-2,  etc-  The  value  of  an  integral 
between  other  integral  limits  can  be  found  by  recalling  .hat 


J~  f(x)dx  = J ^ ( ic ' d x - f f (1  )dx 


Values  of  a . are  eiver  ' Ts.’le  V at  u their  use  con  be  illustrated 
1 


by  an  example. 
Example  6 « 


,1* 


Find  the  £ y dx  given  +h9  table 


xi 

1.5 

1.6 

1 .7 

1 8 

1.9 

2.0 

yi 

48 .09375 

65 . 53600 

87.69705 

115.47360 

149.86915 

192. 
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Solution 


3 

= h £ f(u)du 

= h t £0  « i ’'°-3)  J 


From  Table  V for  six  points  in  the  interval  0-3  we  have 


a . 

i 

.31875 

1.36875 

.71250 

.71250 

- .13125 

.01875 

and  since  h = .1 

i.e 

/ ydx  = .1  (233.721045)  = 23.372104 

1.5 

7.  Solution  of  Differential  Equations 

The  methods  most  frequently  used  in  the  numerical  solution  of 
differential  equations  are  those  called  step-by-step  methods  in  which 
the  values  of  the  dependent  variable  are  calculated  from  a sequence 
of  equally  spaced  values  of  the  independent  variable.  That  is,  if 
we  have  a differential  equation 

(69)  y'  = f(x.y) 

and  the  initial  conditions  x , y , we  compute  a sequence  of  values 

o o 

for  y , say  yx  , y2  , y3  , .....  for  equally  spaced  chosen  values  of 

x , say  Xj^  , Xj  , Xg To  do  this  we  rely  upon  formulas  which 

yield  the  next  value  of  y , say  y^  , from  previously  known  values 

of  y and  its  derivative.  Such  recursion  formulas  for  y^  + 1 b® 

derived  from  the  Lagrangian  formula  for  finding  the  derivative  of  a 
function. 

Let  us  apply  formula  (58)  to  the  case  of  five  points,  we  have 
(See  Table  III) 

i 

* 

| , 

i 

t 
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J2h  ('25yo  + 48y-  * ,i6y-  + 16y3  ‘ ay4  ) 

yi  = iih  ( '3y-'  * i0y-  + 18y= " 6ys  + y‘<  J 


(7°)S  = 


— - — ( y 

iah  v yo 


+ 8y^  - y4  ) 


Bh  < - »0  * "5'-  - + 10V3  + 4».  1 


. K ~ — h ( 3»0  " 1°yi  + ‘ 48,3  + 25,4  > 

If  we  consider  the  foimulas  ; or  y'  , y"  , and  y ' we  can  solve 

I 2 3 

for  y in  terms  oi  these  derivatives  and  y 
J4  •'o 

12h  yx  = -3yQ  - 10y1  + .L«ya  - 6y3  + y^ 


(71) 


12h  y' 


8yl  + 8ya 


a * O ''l  “3  '4 

12h  y'  = - y + 6y,  - 18y,  + lUy  + 3y 

t J3  J0  1 2 3 4 

We  have  from  the  second  equation  of  (71) 

(72)  y^  = yQ  - 8y,  + 8^  - 12h  y' 

Adding  tha  first  and  third  oauation  of  (71)  gives 
12h  ( y + y ) = - 4y  - 4y.  + 4y  + 4y 

-*  Q _ J 4 

y. 


01 


yi  - y3  = 


3h(y'  + y'  ) + y 

W1  3 ' 


and  upon  substitution  of  thx^  into  equation  (72)  we  get 


yA  = y. 


8 [ 


'h(y,'  + y'  ) + y4  ] - 12h  y ' 


Oi 


9y,  = 9y  + 12h(  2y ' - y'  + 2y  ) 
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so  that 

(73)  y4  = y0  * f-  ( 2y;  - va  + 2v;  ). 

This  formula  yields  a value  of  y in  terms  of  the  value  of  y four 
steps  previous  and  the  values  of  the  derivative  of  y at  the  pre- 


4 
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f 


) 

! 

I 

' 


which  may  be  obtained  by  algebraic  manipulation  of  the  Lagrangian 
derivative  formulas  Trial  values  to  start  the  method  of  successive 
approximation  are  given  by  the  simple  Euler  formulas 


(76.) 


+ hv" 
o 


The  procedure  will  be  Illustrated  by  an  example. 


Example  7 . 


Solve  the  differential  equation 


&L  = 

dx 


2x  - y 


Solution  Differentiating  y‘  we  get 


y"  = 2 - y'  = 2 - 2x+y 

and  evaluation  at  (x  ,y  ) yields  y"  = 3. 
The  calculation  may  be  arranged  as  follows 

y'  = 2x  - y 


h/ 24  = .004166667 
yo  = 3. 

hay * / 4 = .0075 
*o 


h/3  = .033333 
4h/3  = .133333 


1.  The  Start  from  Formulas  (75)  and  (76) 


y_i 

3.115 

3.1155 

3.115512 

3-115512 

y' 

-1  3 

-1.315 

-1.3155 

-1.315512 

-i 

yi 

2.915 

2.9145 

2.914513 

2.914513 

-.7 

- .715 

- 7145 

- .714513 

ya 

2.856197 

2.856192 

la ... 

- .456197 

- .456192 

2.  The  Continuation  from  (73)  and  ("4) 


X 

yc 

yc 

yp 

* 

yp 

.9 

3.115512 

-1.315512 

1.0 

3. 

-1. 

1.1 

2.914513 

- .714513 

1.2 

2.856192 

- .456192 

1.3 

2.822456 

-.222456 

2.822463 

- .222463 

1.4 

2.810960 

-.010960 

2.810968 

- .010968 

1.5 

2.019593 

.180407 

2.819600 

.180400 

1.6 

2.846434 

.353566 

2.846440 

.353560 

1.7 

2.889756 

.510244 

2.889763 

-510237 

1.8 

2.947985 

.652015 

2.947991 

.652009 

1.9 

3.019708 

.780292 

3.019715 

.780285 

2.0 

3.103635 

.896365 

3.103641 

.896359 

W-53-35 
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8.  Higher  Derivatives 

The  higher  derivatives  are  obtained  by  repeated  differentiation 

of  equation  (39).  Thus 

(77)  L"  (x)  = L"  (x)  yQ  + Li'  (x)  yx  + + L"  (x)  y n 

where 

(78)  (x)  = [ Pi(,)  ] 

i l dx 

Although  a notation  and  schematic  can  be  devised  for  the  general 
case  it  is  very  unwieldy.  Consequently,  we  shall  develop  only  the 
special  case  of  equally  spaced  intervals  From  equation  (50)  we  have 

(79)  y'(x)  o i 2 L:  (u)y 

n i=o  i i 


so  that 


(80) 

y~(x)  ■ 

= 

n i=o 

L. 

l 

Now 

Li(u)  = 

(-if1 

*T 

n ' / 

“nTT-iTT”  ^U' 

where 

is  the  binomial 

coefficient  notation. 

Let 

(81) 

b.  = 

l 

i '.  (n-i)  - 

then 

L'(u) 

= b. 
i 

d 

du 

n 

[ j— o OPI  J 

and 

(82) 

l;'(u) 

= b. 
i 

_d* 

du 

.n  (u- o ) 

J=0  Jpl 

• 

These  derivatives  may  now  be  evaluated  at  u = 0,  1,  2,  3, 
yield  the  values  of  the  derivatives  at  x = x , xlt  .... 
we  can  write 


(83)  y"  (xk) 


h"Uk)  = -y  UkV„  - *J1y1  + »kV2  - 


where 


. . , n to 

x . Thus 
n 


.+ 


A"  y ] 
kn  n 


hr 


(k) 
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whioh  may  be  compute!  once  and  for  all  and  have  been  tabulated  in 
Table  VI. 

The  continuation  to  higher  dor  natives  is  immediate.  The 
third  derivative  is  given  by 

(84)  y"(xk)  = L~(xk)  = ^(A-;yo  - AJ'1,1  4 ...  4 »~ynJ 

and  the  values  for  k'"  have  been  tabulated  in  Table  VII. 

ki 

Example  8. 

Find  the  first,  second,  and  third  derivatives  at  x = 1 of  the 
function  whose  values  are 


X 

2 

-1 

0 

T 

2 

3 

4 

y 

104 

17 

0 

-1 

8 

69 

272 

Solution 


Using  Tables  III,  VI,  VII  with  n=6,  h = 1 , x3  = 1 , u = 3 
we  write 


y'  (i) 


g-Q-  [ -1(104)  + 9(17)  - 45(0)  + 0(-l)  + 45(8)  - 9(69)  + 1(272)] 
0 • 


y’  (i)  = 


a) 


(2(104)  -2r? (17 1 + 270(01  - 490(-l)  + 270(8)  - 27(69) 

+ 2(272)] 

E- 

= w [3(104)  - 24(17)  + 39(0)  + 0(-l)  - 39(8)  + 24(69) 

* - 3(272)] 


9o  Lagrangian  Multipliers 

Although  the  Lagrangian  multiplier  is  not  directly  associated 
with  the  Lagrangian  formulas  discussed  this  far  it  still  seems 
appropriate  to  include  them  in  this  paper  because  they  are  concerned 
with  the  finding  of  an  extremum  of  a function  which  is  constrained 
by  some  additional  condition.  In  the  principle  of  least  squares 
which  is  so  often  used  in  numerical  analysis  we  are  concerned  with 
the  problem  of  minimizing  the  sum  of  the  squares  of  the  residuals 
of  the  analytical  expression  am  the  values  found  by  substitution  of 
the  known  points.  If  now  an  additional  constraint  is  imposed  upon  tne 
parameters  in  the  analytical  expression  the  problem  could  be  solved 
by  use  of  the  Lagrangian  multiplier. 
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Let  us  consider  a function  of  two  variables;  f(x,y)  and  let  it 
be  desired  to  find  the  maximum  and  minimum  values  of  this  function  if 
furthermore  x and  y are  related  by  g(x,y)  = 0. 

The  necessary  condition  for  an  extremum  is  that  the  total  differ- 
ential vanishes.  Thus  for  the  function  f(x,y) 

(85)  dx  + |-f  Jy  = 0. 

dx  dg 

The  total  differential  of  the  constrained  relation  is 


(86) 


dx  + dy 

3x  3g 


0 • 


If  now  equation  (86) 
X,  and  added  to  equation 


is  multiplied  by  an  undetermined  multiplier, 
(85)  we  have 


<87>  (t§  + xif) dx  + (tv’  + xt§) 


Equation  (87)  would  be  satisfied  if  X 


( 


*1  + \ Is. 

dx  dx 


)- 


(88) 


(8  • >8) 


g(x.y)  = o 


dy  = 0 

were  determined  so  that 


The  multiplier  X is  called  the  Lagrangian  multiplier. 


Let  us  apply  this  multiplier  to  a problem  of  fitting  data  by  the 
principle  of  least  squares.  Consider  the  problem  of  fitting  a set 
of  data  by  the  expression  y = mx  + b . The  principle  of  least 
squares  forms  the  n residuals 

= m + b - (i  = 1,  .....  n) 
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Example  9 

s Fit  a straight  line  to  the  data  such  that  the  slope  is  a minus 
twioe  the  "intercept . 


X 

l 

2 

3 

4 

5 

6 

7 

y 

.4 

1.1 

1.5 

2.2 

2.6 

3.1 

3.8 

Solution 


n = 7,  S1  = 28,  «1  = 14.7,  «2  = 74.1,  S2  = 140 


14.7  - 2(74.1) 
b “ 7 + 4(140)  - 4*^28*)* 

m = .5868 


•133.5 
' 4 5 5 ' 


= -.2934 


y = . 5868x  - .2934 


Che ok : 


X 

l 

2 

_ 3 _j 

4 

5 

6 

7 

y 

.2934 

.8802 

1.467 

2.0538 

2.6406 

3.2274 

3.8142 

I 


This,  of  course,  is  not  the  best  fitting  straight  line;  it  is 
the  best  fitting  straight  line  for  which  m = -2b. 

10.  Accuracy  of  Lagrangian  Formulae 

All  of  the  formulas  derived  in  this  paper  are  based  on  the 
assumption  that  the  given  data  can  be  fitted  by  a polynomial  in  the 
region  under  discussion  and  this  initial  assumption  should  frequently 
be  checked.  The  inherent  error  in  the  Lagrangian  formulae  is  ex- 
hibited in  the  remainder  term  of  equation  (12) 

(n+i ) 

Rn  = P(x) 

where  £ is  an  arbitrary  point  in  the  interval  (0,  n). 

This  remainder  term  may  be  carried  throughout  all  of  the  derivatives 
of  the  formulas  and  is  a measure  of  their  accuracy.  Thus  in  formula 
(58)  for  the  derivative  should  be  added  the  remainder 

hn  (n+i ) 

~ y 

For  the  integrals  we  have  the  expression  (63). 

The  accuracy  should  always  be  tested  whenevor  possible. 
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1922  -24340  271292  -2037253  13329235  -42992730  1022  076  -124133396  1=9123160  399163000 


TABLS  IX. 

Lagrangian  Interpolation  Coefficients 

5 - Points 


n 

A-2 

~A-i 

*0 

A1 

-Aa 

.0000000 

. 0000000 

1 . 0000000 

.0000000 

.0000000 

0.00 

.0008291 

.0065998 

.9998750 

.0067  332 

.0008374 

.01 

.0016493 

.0130654 

.9995000 

.0135986 

.0016827 

.02 

.0024603 

.0193956 

.9983752 

.0205954 

.0025352 

.03 

.04 

.0032614 

.0255898 

.9980006 

.0277222 

.0033946 

.04 

0.05 

.0040523 

.0316469 

. 996d766 

.0349781 

.0042602 

0.05 

.06 

.0048325 

.0376662 

.9956032 

.0423618 

.0051315 

.06 

.07 

.0056016 

.0433468 

.9933810 

.0498722 

.0060079 

.07 

.08 

.0063590 

.0489882 

.9920102 

.0575078 

.0068890 

.08 

.09 

.0071045 

.0544894 

.98 93914 

.0652876 

.00777  40 

.09 

0.10 

.0078375 

.0598600 

.98752  50 

.0731500 

.0086625 

0.10 

.11 

.0085577 

.06506  92 

.9849116 

.0811538 

.0095538 

.11 

.12 

.0092646 

.0701466 

.9820518 

.0892774 

.0104474 

.12 

.13 

.0099580 

.0760814 

.9789464 

.0975196 

.0113425 

.13 

.14 

.0106373 

. 07937  34 

.9755960 

.1058786 

.0122387 

.14 

0.15 

.0113023 

.0845219 

.9720016 

.1143531 

.0131352 

0.15 

.16 

.0119526 

.0890266 

.9681638 

.1229414 

.0140314 

.16 

.17 

.0125879 

.0933870 

.9640838 

.1316420 

.0149266 

.17 

.18 

.0132077 

.0976030 

.9597624 

.1404530 

.0158203 

.18 

.19 

.0138119 

.1016740 

.9652008 

.1493730 

.0167116 

.19 

0.20 

.0144000 

.1056000 

.9504000 

.1584000 

.0176000 

..20 

.21 

.0149716 

.1093806 

.9453612 

.1675324 

.0184847 

.21 

.22 

.0155269 

.1130153 

.9400856 

.1767682 

.0193651 

.22 

.23 

.0160652 

.1165052 

.9346746 

.1651058 

.0202403 

.23 

.24 

.0165862 

.1198490 

.9288294 

.1955430 

.0211098 

.24 

0.25 

.0170898 

.*230469 

.92.  ..8516 

.2050781 

.0219727 

0.25 

.26 

.0175757 

.1260990 

.9166  \ 

.2147090 

.0228233 

.26 

.27 

.0180437 

.1290052 

.9102036 

.2244338 

.0236753 

.27 

.28 

.0184934 

.1317658 

.9035366 

.2342502 

.0245146 

.28 

.29 

.0189248 

.1343806 

.8966432 

.2441564 

.0253437 

.29 

0.30 

.01933-5 

.1368500 

.88952  50 

.2541500 

.0261625 

0.30 

.31 

.0197314 

.1391740 

.8821838 

.2642290 

.0269701 

.31 

.32 

.0201062 

.1413530 

.8746214 

.2743910 

.0277658 

.32 

.33 

.0204619 

.1433870 

.8668398 

.2846340 

.0285436 

.33 

.34 

.0207981 

.1452766 

.8588408 

.2949554 

.0293179 

.34 

0.35 

.0211148 

.1470219 

.8506266 

.3053531 

. 0300727 

.35 

.36 

.0214118 

.1486234 

.8421990 

.3158246 

.0308122 

.36 

.37 

.0216890 

.1500814 

.8335604 

.3263676 

.0315355 

.37 

.38 

.0219461 

.1513966 

.8247128 

.3369794 

.0322419 

.38 

.39 

.0221832 

.1525692 

.8156586 

.3476578 

.0329303 

.39 

0.40 

.0224000 

.1536000 

.8064000 

.3584000 

.0336000 

0.40 

.41 

.0225965 

.1544894 

.7969394 

.3692036 

.0342500 

.41 

.42 

.0227725 

.1552382 

.78727  92 

.3800658 

.0348795 

.42 

.43 

.0229281 

.1558468 

.7774220 

.3909842 

.0354874 

.43 

.44 

.0230630 

.1563162 

.7673702 

.4019558 

.0360730 

.44 

0.45 

.0231773 

.1566469 

.7571266 

.4129781 

.0366352 

.45 

.46 

.0232709 

.1568398 

.7466936 

.4240462 

.0371731 

.46 

.47 

.0233438 

.1568956 

.73607  42 

.4351634 

.0376857 

.47 

.48 

.0233958 

.1568154 

.7252710 

.4463206 

.0381722 

.48 

.49 

.0234271 

.1565998 

.7142870 

.4575172 

.0386314 

.49 

0.50 

.0234375 

.1562500 

.7031250 

.4687500 

.0390625 

0.50 

-P 

Aa 

H 

< 

1 

Ao 

A-1 

” A-  2 

\m 

TABU  U. 


LagranKian  Interpolation  Coefficients 

5 - Points 


B 

m 

Ao 

| 

■9 

m 

0.50 

.0234375 

.1562500 

.7CK51Iit>0 

.4667500 

.0390625 

0.50 

.51 

.0214271 

.1557..68 

.6917  inO 

.4800162 
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0.55 
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.6202650 
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.57 

.58 

.0227725 

.1487422 

.6077912 
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.0413755 

.58 

.59 

.0225965 
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.0415070 

.59 
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.1456000 
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.0416000 
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.0221832 

.1438462 
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.62 
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.1419726 
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.1399804 
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.63 

.64 
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.1378714 
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.0211148 

.1356469 

.5165016 

.6394781 

.0414477 

0.65 
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.0207981 

.1333086 

.6029368 
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.0412859 

.66 

.67 

.0204619 

.1308580 

.4892528 

.6622210 

.0410776 

.67 

.68 
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.1282970 

. 4754534 

.6735590 

.0408218 

.68 

.69 

.0197314 

.1256270 

.4615428 

.6848700 

.0405171 

.69 

0.70 

.0193375 

.1228500 

.4475250 

.6961500 

.0401c>25 

0.70 

.71 

.0189248 

.1199676 

.4334042 

.7073954 

.0397567 

.71 

.72 

.0184934 

.1169818 

.4191846 

.7 166022 

.0392936 

.72 

.73 

.0180437 

. 113*942 

.4048706 

.7297668 

.0387868 

.73 

.74 

.0175757 

.1107070 

.3904604 

.7408850 

.0382203 

.74 

0.75 

.0170898 

.1074219 

.3759766 

.7519531 

.0375977 

0.75 

.76 

.0165062 

.1040410 

.3614054 
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. 76 

.77 

.0160652 
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.0155269 
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.0144000 
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.3024000 
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0.80 
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.2874918 
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.82 

.0132077 
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.2725304 
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.8382060 
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.0266867 
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.0252915 
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.3894054 
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.8993348 

.0222808 

.89 
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.0078375 
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.0206625 

0.90 

.91 

.0071045 

.0432964 

.1363124 

.9188466 

.0189670 

.91 

.92 

.0063590 

.0386842 

.1210982 

.9284198 

.0171930 
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.93 

.0056016 

.0340158 

.1058880 
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.0153389 

.93 
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.0048325 

.0292942 
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.9471778 

.0134035 

.94 
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.0040523 

.0245219 

.0755016 

.9563531 

.0113852 

0.95 

.96 

.0032614 

.0197018 

.0603366 
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.0092826 
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.0148366 
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.0070942 
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.0016493 

.0099294 

.0300920 
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.0049828 
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Differentiation  Coefficients  - Fractional  Form 
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a' 
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A' 
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a/ 
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\ 

ko 
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k2 

ki 

k« 
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kc 

k7 

kp 
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n = 2 (3  points)  Dn  - 2 


0 

-3 

4 

-1 

! 2 

1 

-1 

0 

1 

LU 

n = 3 (4  points)  On  = 6 


0 

-li 

13 

-9 

2 

3 

I 

_ o 

-3 

6 

-1 

2 

n = 4 (5  points)  Dn  = 12 


— r— n 

0 

43 

-36 

16 

-3 

4 

1 | -3 

-10 

18 

-6 

1 

3 

— 

-8 

0 

8 

-1 



n = 5 (6  points)  Dn  = 60 


0 

-137 

3C0 

-300 
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-75 

1? 

5 

1 

-12 

-6  5 
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-60 

20 

_ *T 

4 

2 

3 

-30 

-20 

•'0 

-15 

2 

3 

n = 

S (7  points)  Dn 
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0 
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-225 

72 
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6 

1 
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-77 
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50 
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2 

5 

2 

2 

-24 

-35 

’0 

-30 

i 

-1 

4 

3 

-1 

9 
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0 

45 

-9 

1 

n = 7 (8  points)  Dn  = 420 


0 
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60 
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1 
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-609 
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54 
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n 

c 

10 
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-350 

140 
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4 

5 

T 

-4 

42 

-252 

-105 
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2a 

-3 

n = 8 (9  points)  Dn  = 840 


0 

1 

2 

3 

4 

-2283 
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15 

-5 

3 
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-1338 
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60 

-32 
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2940 

-793 

-420 
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-2940 

1680 

-378 
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-14700 

2840 
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0 

9104 

-1470 
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-420 
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48 
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32 

— 
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15 

-5 

3 

-3 

d 

7 

6 

5 
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n = 9 

(10  points)  Dn 
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0 

-7129 
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-45360 

70560 

-79300 
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-35280 

12960 

-2835 
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9 

1 

-280 

-4329 

10080 

-11760 
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-8820 

4704 

-1680 
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-35 

8 

2 

35 
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-2754 

5030 

-4410 

2940 

-1470 
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10 

7 

3 

-10 
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-1554 

3780 

-1390 
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54 

-5 

6 

4 

P 
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240 

-45 

4 

5 

k 
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-K 
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A, 
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k 
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tabu:  iv. 


Differentiation  Coefficients  - Decimal  Form 


! 

0 -1. -33333 

3.000000 

-1 

.500000 

.333333 

» j -.  33333 

- . 500000 

1 

. 000000 

- . 166667 

-2.063333 

4 . 000000 

-3 .000000 

1.333333 

- .250000 

- .250000 

-.833333 

1 . 500000 

- . 500000 

.0*3333 

.033333 



- .666667 

0 

.666667 

- .083333 

0 

-2.283333 

5 . 000000 

-5.000000 

3.333333 

— 

-1.250000 

.200000 

1 

- . 200000 

-1 .0“3333 

2.000000 

-1.000000 

.333333 

-.050000 

2 

.050000 

- . 600000 

-.333333 

1.000000 

-.210000 

.C  33333 

0 

-2.450000 

6.000000 

-7 .500000 

6 .666667 

-3.760000 

1.200000 

- .166667 

1 

-.156667 

-1 .283333 

2.500000 

-1.666667 

.833333 

- .250000 

.033333 

2 

.033333 

-.400000 

- . 5^3333 

1.333333 

- .500000 

.133333 

- .016667 

3 

-.016667 

.150000 

- .7  50000 

0 

.750000 

-.150000 

.016667 

0 

-2.5928  57 

7 . 000000 

-10.500000 

11.666667 

-8  .7  50000 

4.200000 

-1  .166667 

.1428  57 

1 

- .142857 

-1.450000 

3.000000 

-2.500000 

1 . 666667 

- .750000 

.200000 

-.023810 

2 

.023810 

-.333333 

- .783333 

1.666667 

- .833333 

.333333 

-.083333 

.009524 

3 

- .009524 

. 100000 

- .600000 

-.250000 

1 .000000 

-.100000 

.066667 

-.007143 

I 


-2.717857 

8.000000 

-14.000000 

18 .666667 

-17.500000 

11.200000 

-4.666667 

1.142857 

-.125000 

-.125000 

-1.592857 

3.500000 

-3 . 500000 

2.915667 

-1 .750000 

.700000 

-.166667 

.017857 

.017857 

- .285714 

- .950000 

2.000000 

-1.250000 

.666667 

- .250000 

.057143 

-.005952 

-.005952 

.071429 

-.500000 

-.450000 

1.250000 

- . 500000 

.166667 

-.035714 

.003571 

.003571 

-.038095 

.200000 

-.800000 

0 

. 800000 

- .200000 

.038095 

-.003571 

-2.328968 

9.000000 

-18.000000 

28 .000000 

-31.500000 

25.200000 

-.111111 

-1.717857 

4.000000 

-4.666667 

4.666667 

-3.500000 

.013889 

-.250000 

-1.092857 

2.333333 

-1.750000 

1 .166667 

-.003968 

.053571 

- .428571 

-.616667 

1 . 500000 

- .750000 

.001984 

-.023810 

.142857 

-.666667 

- .200000 

1 . 000000 

1 .666667 
- .583333 
.333333 


5.142857 

-1.125000 

.111111 

-.666667 

.142857 

-.013889 

.200000 

-.041667 

.003  968 

-.107143 

.021429 

-.001984 

.095238 

1 

-.017867 

.001587 

f 

' 
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TABU’  V. 

Integration  Coefficients  ) 


n = 2 (3  Points) 


n = 4 (5  Points) 


EBB 

0 - 1 

0 - 2 

0 - 3 

0-4 

.348611111 

.322222222 

.337600000 

.311111111 

M 

.897222222 

1.37  7777778 

1.275000000 

1.422222222 

- . 366666667 

.266666667 

. 900000000 

.633333333 

ESli 

.147222222 

.044444444 

.525000000 

1.422222222 

- .026368889 

- .011111111 

-.037500000 

.311111111 

n = 5 

( 6 Points) 

Hi 

1 

o 

— 

CVJ 

l 

o 

0 - 3 

0 - 4 

0 - 5 

1 

.329861111 

.990972222 

- .554166667 
.334722222 

- .120133889 
.0187  50000 

.311111111 
1.433333333 
.156555556 
.156555560 
- . 066666667 
.011111111 

.319750000 
1.363750000 
.712600000 
.712500000 
- .131250000 
.018750000 

.311111111 
1 .422222222 
.533333333 
1.422222222 
.311111111 
0 

.329861111 

1.302083333 

.068055556 

.368055556 

1.302083333 

.329861111 

n = 6 (7  Points  ) 


EH 

0 - 1 

CQ 

1 

o 

0-3 

0-4 

0 - 5 

0 - 6 

ao 

“3 

as 

«c 

.315591931 

1.076587301 

- .768204365 
.620105820 

- .334176587 
.104365079 

-.014269179 

.301322751 
1.492063492 
.003730169 
.351322751 
- .213492064 
.069841270 
- .009788359 

.305903571 
1.446428571 
.518303571 
.971428571 
- .325446428 
.096428571 
- .012946428 

.302645503 
1.473015873 
.406349206 
1 .591534392 
.184126984 
.050793651 
- .008465609 

.307126323 
1 . 439492C63 
.527033730 
1.322751323 
.961061508 
.466269841 
- .022734788 

.292857143 

1.542857143 

.192857143 

1.942857143 

.192857143 

1.542357143 

.292857143 

n = 7 ( 8 Points  ) 


0 - 1 

0-2 

ISKHEBi 

0 

1 

HK3EKH 

.304^24537 

.292857143 

.295758928 

.294179894 

.295758928 

.292357143 

.304224537 

1.156159060 

1.551322751 

1.516741071 

1 .532275132 

1.513063822 

1.542857143 

1.449016203 

-1.006919642 

-.169047619 

.307366071 

.228571428 

.288318452 

.192857143 

.535937500 

1.017964616 

.647619047 

1.322991071 

1 .937330687 

1.720610119 

1.942357143 

1.210821759 

- .732035383 

-.509738359 

- .677003928 

- .112169312 

.563202711 

.192857143 

1.210821759 

.343080357 

.247619048 

.307366071 

.228571428 

.704985119 

1.542857143 

.535937500 

- .093840939 

-.069047619 

- .03325d92d 

- .067724807 

- .102306547 

.292857143 

1.449016203 

.011367394 

.003465608 

.010044643 

.008465608 

.011367394 

0 

.304224537 
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TABUt  VI. 

Second  Derivative  Coefficients 


W-53-35 
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